{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "configured-penetration",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.07987877, 0.04438157, 0.08815771, 0.06126131, 0.1400486 ,\n",
       "       0.08051048, 0.05087864, 0.03734621, 0.14670686, 0.13679783,\n",
       "       0.11039573, 0.07017681, 0.0429424 , 0.37910509, 0.15169127,\n",
       "       0.08801036, 0.10182255, 0.09908553, 0.11638002, 0.09067552,\n",
       "       0.08478104, 0.09664813, 0.07846983, 0.17283006, 0.10483507,\n",
       "       0.17483697, 0.10958365, 0.06700379, 0.10410579, 0.04921869,\n",
       "       0.14325543, 0.09816408, 0.08652314, 0.06651088, 0.0464032 ,\n",
       "       0.07998701, 0.06426596, 0.04832661, 0.10648687, 0.07773626,\n",
       "       0.05121746, 0.16022885, 0.08096374, 0.08985915, 0.05227419,\n",
       "       0.04702343, 0.09171804, 0.09419053, 0.16234071, 0.1719272 ,\n",
       "       0.0533473 , 0.1038362 , 0.15054845, 0.2134348 , 0.09346092,\n",
       "       0.14740421, 0.09349584, 0.13446947, 0.12710528, 0.14236995])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pymc3 as pm\n",
    "\n",
    "from data_generation import Long_Data_gener\n",
    "import numpy as np\n",
    "from nonparameter_classfication import NonparaClassfication\n",
    "import parameter_classfication\n",
    "import pandas as pd\n",
    "from tqdm import tqdm\n",
    "import scipy.special as special\n",
    "import theano.tensor as tt\n",
    "\n",
    "\n",
    "\n",
    "t=3\n",
    "k=3\n",
    "i=20\n",
    "N=200\n",
    "# 实例化对象\n",
    "gener = Long_Data_gener()\n",
    "classf = NonparaClassfication()\n",
    "# 生成被试高阶能力\n",
    "theta = gener.m_theta(t, mean=0, growth_rate=0.3, n=N, cor=0.6)\n",
    "\n",
    "# 生成被试属性\n",
    "beta_theta = np.random.normal(1,0.25,size=(k,t))\n",
    "_lambda__ =np.random.normal(-0.5,0.25,size=(k,t))\n",
    "att, z = gener.attribute(theta, beta_theta,_lambda__, return_possi=True)\n",
    "# 生成题目Q矩阵\n",
    "Q = gener.Q_generate_long(t, i, k)\n",
    "# 生成模拟作答\n",
    "cov = np.array([[0.35,-0.35*0.5],[-0.35*0.5,0.35]])\n",
    "item_para = np.random.multivariate_normal([4.394449154672438,-2.197224577336219],cov=cov,size=(i*t))\n",
    "lambda_k = item_para[:,0]\n",
    "lambda_0 = item_para[:,1]\n",
    "resp = gener.simu_resp(t, att, Q, lambda_k=lambda_k, lambda_0=lambda_0)\n",
    "\n",
    "\n",
    "\"\"\"被试数据生成到此结束\"\"\"\n",
    "# 获取被试所有可能的属性掌握模式\n",
    "all_pattern =gener.all_pattern(k)\n",
    "# 计算所有可能的属性掌握模式对应的理想作答\n",
    "class_ideal_resp = gener.long_dina_ideal_resp(t,Q,all_pattern)\n",
    "\n",
    "np.exp(lambda_0)/(1+np.exp(lambda_0))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "vietnamese-manor",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.394449154672438"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "item_para.shape\n",
    "special.logit(0.9)-special.logit(0.1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "jewish-agenda",
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as long_dina:\n",
    "    theta_mu = pm.Normal(\"theta_mu\",0,1,shape=t-1)\n",
    "    x = tt.zeros((t))\n",
    "    theta_mu = tt.inc_subtensor(x[1:],theta_mu)\n",
    "    _ = tt.zeros((t,t))\n",
    "    inv_diag_cov = pm.Gamma(\"inv_diag_cov\",1,1,shape=t-1)\n",
    "    inv_dtriangel_cov = pm.Normal(\"inv_dtriangel_cov\",0,1,shape=int((t**2-t)/2))\n",
    "    _ = tt.inc_subtensor(_[0,0],1)\n",
    "    cnt = 0\n",
    "    for t_ in range(1,t):\n",
    "        _ = tt.inc_subtensor(_[t_,t_],inv_diag_cov[t_-1])\n",
    "        for t__ in range(t_-1):\n",
    "            _ = tt.inc_subtensor(_[t_,t__],inv_dtriangel_cov[cnt])\n",
    "            cnt+=1\n",
    "    y = pm.math.matrix_inverse(_.dot(_.T))\n",
    "    theta_cov_mat = pm.Deterministic(\"theta_cov_mat\", y)\n",
    "    theta_ = pm.MvNormal(\"theta\",theta_mu,theta_cov_mat,shape=(N,t))\n",
    "    beta = pm.Lognormal(\"beta\",0,0.25,shape=(k,t))\n",
    "    lambda_ = pm.Normal(\"lambda_\",0,0.25,shape=(k,t))\n",
    "    Bound_normal = pm.Bound(pm.Normal,0)\n",
    "    beta_item = pm.Normal(\"beta_item\",-1.096, 0.25,shape=(t,i))\n",
    "    delta_item = Bound_normal(\"delta_item\",0,0.25,shape=(t,i))\n",
    "    att_ = tt.zeros((t,N,k))\n",
    "    for t_ in range(t):\n",
    "            x = tt.outer(theta_[:,t_],beta[:,t_])+lambda_[:,t_]\n",
    "            x = tt.exp(x)/(1+tt.exp(x))\n",
    "            att_ = tt.inc_subtensor(att_[t_],x)\n",
    "    att_ = pm.Deterministic(\"att_possi\",att_)\n",
    "    att_student = pm.Bernoulli(\"att_student\",att_,shape=(t,N,k))\n",
    "\n",
    "    resp_ = tt.zeros((t, N, i))\n",
    "    for t_ in range(t):\n",
    "        for i_ in range(i):\n",
    "            x = att_student[t_]>=Q[t_,i_]\n",
    "            eta_ = tt.prod(tt.cast(x, dtype=\"int32\"), axis=1)\n",
    "            resp_ = tt.inc_subtensor(resp_[t_,:,i_],tt.exp(beta_item[t_,i_]*eta_+delta_item[t_,i_])/(1+tt.exp((beta_item[t_,i_]*eta_+delta_item[t_,i_]))))\n",
    "    pm.Bernoulli(\"response\",resp_,observed=resp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "rational-combination",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (2 chains in 4 jobs)\n",
      "CompoundStep\n",
      ">NUTS: [delta_item, beta_item, lambda_, beta, theta, inv_dtriangel_cov, inv_diag_cov, theta_mu]\n",
      ">BinaryGibbsMetropolis: [att_student]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='8007' class='' max='10000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      80.07% [8007/10000 4:26:30<1:06:20 Sampling 2 chains, 562 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# step1 = pm.NUTS(target_accept=0.95)\n",
    "\n",
    "with long_dina:\n",
    "#     step = [pm.NUTS([delta_item, beta_item, lambda_, beta, theta_, chol, theta_mu], target_accept=0.95),\n",
    "#        pm.BinaryMetropolis([att_student],transit_p=0.9)]\n",
    "    trace = pm.sample(3000,tune=2000,chains=2,cores=4)\n",
    "# >NUTS: [g, s, lambda_, beta, theta, theta_cov_LJK, theta_mu]\n",
    "# >BinaryGibbsMetropolis: [att_student]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "careful-slovak",
   "metadata": {},
   "outputs": [],
   "source": [
    "import arviz as az"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "grave-illness",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = az.summary(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "greater-content",
   "metadata": {},
   "outputs": [],
   "source": [
    "az.plot_energy(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "electoral-tactics",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"ACCR:\",np.mean((trace.get_values(\"att_student\").mean(axis=0)>0.5)==att))\n",
    "print(\"PCCR:\",np.mean(np.prod((trace.get_values(\"att_student\").mean(axis=0)>0.5)==att,axis=2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "improved-palmer",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.mean(np.product(np.product((trace.get_values(\"att_student\").mean(axis=0)>0.5)==att,axis=0),axis=1))\n",
    "print(\"Long_PCCR:\",x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "overhead-episode",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.mean(az.summary(trace,\"beta_item\")[\"mean\"]-lambda_k))\n",
    "print(np.mean(az.summary(trace,\"delta_item\")[\"mean\"]-lambda_0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "european-remains",
   "metadata": {},
   "outputs": [],
   "source": [
    "az.summary(trace,\"inv_dtriangel_cov\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "invisible-african",
   "metadata": {},
   "outputs": [],
   "source": [
    "az.summary(trace,\"theta_cov_mat\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "painful-establishment",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(trace.get_values(\"theta\").mean(axis=0)-theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "continuing-cabinet",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "criminal-theta",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext watermark\n",
    "%watermark -n -u -v -iv -w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adapted-junior",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
